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We use numerically exact quantum Monte Carlo (QMC) to compute the properties of three 
dimensional bipolarons for interaction strengths where perturbation theory fails. For intermediate 
electron-phonon coupling and Hubbard U, we find that bipolarons can be both small and light, a 
prerequisite for bipolaron superconductivity. We use the QMC results to make estimates of transition 
temperatures, which peak at between 90— 120 K and are demonstrated to be insensitive to Coulomb 
repulsion and impurities. 



I. INTRODUCTION 

The mechanisms of three-dimensional superconduc- 
tors with high transition temperatures (such as the bis- 
muthates) are widely acknowledged to have their ori- 
gins in the electron-phonon interaction. In at least some 
materials, a proposed mechanism involves the pairing 
of phonon dressed electrons (polarons) to form bipo- 
larons. At sufficiently low temperatures, these bipolarons 
form a Bose-Einstein condensate with superconducting 
properties^. The main criteria for bipolaronic super- 
conductivity with significant transition temperatures are 
bipolaron mobility (low effective massjP with sufficient 
density (small bipolarons). Although small (bi)polarons 
are often viewed as immobile states that may be easily 
localized by disorder, small mobile bipolarons that can 
lead to superconductivity with high transition tempera- 
tures have been found in 1 and 2 dimensions!^' 

Polarons occur naturally in nearly all media, from plas- 
mas and ultra cold atoms through normal bulk mat erials 
and possibly high-temperature superconductors!^ Lan- 
dau introduced the polaron concept in 1933 to study lat- 
tice polarization due to the motion of electrons through 
ionic solids^ An electron moving though a crystal lattice 
creates distortions that follow its trajectory, producing a 
phonon cloud that propagates though the system, sur- 
rounding the electronPEl 

Bipolarons are formed when two polarons interact with 
each other using phonon-mediated interactions to form 
pairs. These pairs can be strongly bound and travel 
though the lattice as a single composite particle. There 
are many different interaction types that lead to bipo- 
laron creationP In this paper, we consider an extended 
Hubbard-Holstein model, where the Coulomb interaction 
is purely on-site, and an extended-Holstein interaction 
couples the electron density to lattice vibrations on the 
same site and near-neighbor sites, simila r to the inter- 
action introduced by Bonca and TrugmanP^l Increasing 
this inter-site interaction allows electrons to form bipo- 
larons and overcomes the Coulomb repulsion so that large 
bipolarons become local (small) bipolarons with approx- 
imately the size of the lattice constant.^ The site- local 
Holstein model describes an extreme short-range limit 
where electrons can form pairs on single atoms. In this 
case, the electron-phonon interaction has to overcome the 



on-site repulsion to form a bipolaronP^ We note that the 
interactions used here are distinct from the specific forms 
required to describe bismuthate superconductors.^ 

Recent studies of low-dimensional bipolarons have uti- 
lized various numerical methods and analytical tech- 
niques. One-dimensional bipolarons have been found 
to be relevant in describing strong electron-phonon in- 
teractions in low-dimensional organic semiconductors,^ 
and it is possible that three-dimensional (3D) bipo- 
larons in a strong magnetic field simplify into one- 
dimensional (ID) bipolarons. 14 Two-dimensional (2D) 
bipolarons have been investigated extensively in the 
study of two dimensional conductors and high temper- 
ature superconductors such as the cupratesP^ 

Properties of the short range Hubbard-Holstein bipo- 
laron model have been established on small lattices us- 
ing exact diagonalizatiorP^ and an optimized approach 
for exact diagonalization at weak coupling! 1 -^ Advanced 
variational techniques,^ density-matrix renormalization- 
group,^ and various quantum Monte Catlo (QMC) 
approaches^EH have all been used to study 2D bipolaron 
systems. It is found in one dimension that only a small 
attractive force between electrons leads to pairing, the 
most important factor being the nearest neighbo r int erac- 
tion when considering long ranged interactions ! 7 * 19 ! Two 
dimensional work also shows us that inclusion of nearest- 
neighbor interaction is responsible for significant change 
compared with the Hubbard-Holstein model. 8 20 Light 
bipolarons are found on a simple square lattice, showing 
that elaborate lattices are not needed to create small light 
pairs that have the potential to form Bose condensatesP 
On change of dimension from ID to 2D, an increase in 
electron-phonon coupling constant and nearest neighbor 
attraction is needed to create on-site bipolaronsP^ 

There have been several notable studies of 3D bipo- 
larons, especially in relation to superconductors, to un- 
derstand the reasons why cuprates and oth er lay ered 
superconductors are different to 3D materials! 21 * 22 * Sev- 
eral publications have examined the differences between 
the binding of bipolarons in two and three dimen- 
sions, concluding that a greater attra ction is needed to 
form stable bipolarons on 3D latticesPMl It has been 
found that the probability of bipolaron formation in- 
creases with decreasing dimensions or increase in the 
crystal anisotropy.^ Variational studies of the region 
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of existence of the three-dimensional singlet bipolaron 
have allowed investigation of the relationship between 
the critical value of the electron-phonon coupling con- 
stant and the dielectric properties of the medium,^ 
concluding that conditions in alkali-halide crystals are 
not suitable even for metastable bipolarons and that 
three-dimensional continuum bipolarons do not exist in 
La2Cu04. However, metal-ammonia systems potentially 
lie within the region of existence for three-dimensional 
continuum bipolaronsP^ In the context of 3D polarons, 
we have also examined binding to attractive impuri- 
ties, showing that polarons are localized when the im- 
purity potential is around four effective hoppings in 
magnitude.^ 

The work presented here goes beyond previous work 
by considering exact solutions for 3D extended Hubbard- 
Holstein bipolarons (both numerical and analytic). Ex- 
act solutions arc important to understand regions of 
the parameter space where perturbative approximations 
break down, and as we will show are essential for under- 
standing the regions where bipolaronic superconductivity 
is strongest. The paper is organized as follows: In Sec. 
PH we introduce the model. The Lang-Firsov transfor- 
mation is performed before a brief overview of the con- 
tinuous time quantum Monte Carlo simulation method. 
Sec. |III| presents quantum Monte Carlo results for sin- 
glet bipolaron properties, including total energy, number 
of associated phonons, inverse mass, and average bipo- 
laron size. Finally in Sec IV we look at the possibil- 



ity of Bose Einstein condensation of three-dimensional 
bipolarons. For completeness, in the Appendix, we con- 
sider the U — V model corresponding to the high phonon- 
frequency limit. Solving the equation analytically we find 
the binding conditions for varying on-site attraction U 
and nearest neighbor repulsion V . 



II. MODEL AND METHODS 
A. Model 

The extended Hubbard-Holstein model used here has 
its basis in a general electron-phonon Hamiltonian with 
elec tros tatic repulsion, which is written in the following 
fornPl 



The first term in the equation expresses the kinetic en- 
ergy of electrons moving from site to site. The element 
t is the hopping integral for an electron moving between 
neighboring sites. 

The second term in the Hamiltonian represents the 
Coulomb repulsion v between two electrons. Here the re- 
pulsion is approximated to have the Hubbard form, and 
long ranged interactions are assumed to be insignificant 
due to screening in the material so the repulsive term 
has the form flHubbard = n »t"q, where U is the 

magnitude of the repulsiorPET Note that near-neighbor 
interactions are in principle allowed, but neglected in sec- 
tion |III| and appendix [3] of this paper. 

The final three terms include the effects of lattice vi- 
bration. The ion momentum is described by the P m op- 
erator and the ion displacement is signified by Here 
we take t; m to be one dimensional, which is an approxi- 
mation that could relate to phonon modes polarized by a 
strong electric field, a three-dimensional molecular crys- 
tal with molecular ordering along a single direction or 
possibly radial phonon modes. Following Ref. 8j we take 
the force function to be, 



(2) 



This describes interaction between electrons on sites at 
vectors n and vibrating ions between valance sites at 
positions m. li are the vectors between nearest neighbor 
valence sites at r and r' . An effective electron-electron 
interaction can be defined as 



*Ar[r,r'] =J2f m [r]f m+Ar [r'] 



(3) 



so that for the chosen force function, &o[r, r']/$ [0, 0] = 
7 where the nearest neighbor interaction strength (7) 
strictly has the value 1/z. The reason for Ar will be 
explained later on in the paper. Here we will also modify 
7 to investigate the effects of turning on the inter-site 
interaction. For large phonon frequency, this interaction 
can then be mapped directly onto a U — V model. For 
7 = 0, a Holstein interaction is recovered, equivalent to 
fm{ n ) — K fimn- The shift in 7 is equivalent to moving 
the vibrating ions within the unit cell so that they get 
closer to the sites host electrons. 



H = -t E C n> ,a C ™,<> 
< rin 7 >,cr 



.) ) c na c na c ii' a ' C n'cr' ■ / , o y 

nn'acr' m 



pi 

E rn 



(i) 



where n and m represent vectors to electrons and ions 
respectively, c (cH are the creation (annihilation) oper- 
ators for electrons, M is the ion mass, u> is the phonon 
frequency and a is the z component of the electron spin. 



B. Lang-Firsov transformation 

In the limit that phonon frequency becomes infinite; 
the model described in Eq. ([I]) can be mapped onto a 
U — V model, consisting of an on-site Hubbard U and 
inter-site Hubbard V. The mapping uses a Lang-Firsov 
canonical transformation^, which creates a new Hamil- 
tonian H = er s He s and wavefunction = e |"0) ^ 
where H = H + [S, H] + [S, [S, #]] + ••• and S = 
gn(d' — d). Here g is a dimensionless interaction con- 
stant proportional to the force and d'(d) is the phonon 
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creation (annihilation) operator. Under this transforma- 
tion, the creation operators for electrons and phonons 
become, 



c\ ->• c\ = cjexp 



3 

dt^rft =d t + ^g.. n . 



On transforming the atomic Hamiltonian (t - 
electron and phonon subsystems are decoupled: 



(4) 
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H at = - £ W E III + ^ E + i),(8) 



The function $ and a dimensionless interaction param- 
eter A = Ep/W are introduced to simplify the Hamil- 
tonian, where W is the half band-width zt, and E p = 
J2jfoj/ 2Mu)2 = $o(0,0)/2Mw 2 is the polaron shift, 
leading to: 
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Transformation of the tight-binding Hamiltonian leads to 
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Where = /,j /lo\/2Muj. When the phonon frequency 
is very large, the ground state contains no real phonons, 
and there is a further simplification leading to a modified 
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act solutions of the transformed Hamiltonian in the large 
phonon frequency limit can be found in the Appendix for 
comparison with the numerical results. 



C. Computational methods 

We use the continuous-time quantum Monte Carlo 
method, which has been used to simulate the screened 
Hubbard-Frohlich bipolaron in ID and 2E02S1. A more 
in-depth overview of our algorithm has been presented 
in a previous paper and so will not be repeated^. 
The continuous-time quantum Monte Carlo algorithm is 



based on path integrals, where each path rj(r) exists in 
imaginary time and represents a single particle in the 
system. The algorithm probes path configurations which 
are each assigned a weight exp(A) where: 
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Here Ar = r(j3) — r (0) is the distance between the end 
points of the paths in the non-exchange configuration, 
the phonon frequency Co — hu/t, and inverse temperature 
f3 = t/ksT. i = 1, 2 and j = 1, 2 represent the fermion 
paths. v(r\ 1 r2) = US ri , r2 is an instantaneous Hubbard 
repulsion between electrons. The paths lie between the 
range r € 0, /3, and are formed from straight segments 
punctuated with 'kinks' representing site to site hopping. 
The algorithm is used to compute (bi)polaron energy, the 
number of phonons in the system, the effective mass of 
the (bi)polaron and the radius of the bipolaron. 



III. QUANTUM MONTE-CARLO AND 
INTERMEDIATE PHONON FREQUENCY 

Using a continuous-time Quantum Monte Carlo code 
we simulated the extended Hubbard-Holstein model with 
nearest neighbor interaction strengths of 7 = 0, 0.25 and 
0.5 on a cubic lattice. The magnitude of the nearest- 
neighbor component of the electron-phonon interaction 
has been shown in both one and two dimensions to 
be the m ost significant contributing factor to bipolaron 
properties^'. The simulation produces exact numerical 
solutions for the total energy, average number of excited 
phonons, mass and size of singlet bipolarons. 

In the following we examine only the singlet bipo- 
laron for a range of U/t and A at inverse temperature 
/3 = 14 with fixed hu/W = 1, where the half band width 
W = 6t, which is towards the lower end of the interme- 
diate phonon frequency limit. Calculations are carried 
out on an infinite lattice where particles are confined to 
within 50 lattice spacings of each other, which is suffi- 
ciently large to avoid the majority of finite size effects. 
The most sensitive property to finite size effects is the 
bipolaron radius, which does not become infinite when 
the polarons are not bound into a bipolaron. All errors 
are determined using bootstrap re-sampling on the sim- 
ulation data and are displayed as 3 standard errors. 
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FIG. 1. Total ground state energy of a singlet Hubbard- 
Holstein bipolaron simulated by CTQMC (panel a) and sin- 
glet bipolaron with nearest-neighbor interaction strength of 
7 = 0.25 (panel b) and 7 = 0.5 (panel c). Similar to large 
phonon frequency, it takes significant negative Hubbard U to 
bind on-site bipolarons. For large inter-site interaction and 
large A, a crossover is seen between on-site and inter-site bipo- 
larons. 



The binding of bipolarons can be established from the 
total energy. Figure [l|a) depicts the total energy cal- 
culated from our Monte Carlo calculations when 7 = 



(Holstein interaction). Diagonal lines show the presence 
of pairs of electrons on a single site. As the Hubbard on- 
site repulsion U is increased, the on-site pairs are pushed 
apart creating pairs of polarons. As the electron-phonon 
coupling constant A is increased we see that larger Hub- 
bard U is needed to break apart the on-site pairs. When 
the pair is unbound the energy of the bipolaron does not 
change with respect to increasing U, and the line is hori- 
zontal. A key difference here is that it takes a large neg- 
ative Hubbard U to bind on-site bipolarons in contrast 
to ID and 2D systems. 

Plots of the total energy of the bipolaron formed 
when the electron-phonon interaction contains a nearest- 
neighbor component of 7 = 0.25 and 7 = 0.5 are also 
shown [Figs, [ljb) and [Tjjc) respectively]. It is imme- 
diately apparent that the rapid transition from on-site 
bipolaron to free polarons is smoothed out with the ad- 
dition of nearest-neighbor interaction. There is no dra- 
matic change between the total energy of the Holstein 
and extended Holstein bipolarons for electron-phonon 
coupling A = 0.2, as there is insufficient inter-site inter- 
action to bind an off-site bipolaron. This is in contrast 
to ID and 2D where significant qualitative changes to all 
bipolaron properties are found when inter-site interaction 
is switched on. 

We observe that the U value corresponding to the point 
of inflection is reduced with increasing 7. This does not 
correspond to bipolarons which are more weakly bound 
(i.e. easier to unbind when the Hubbard U is switched 
on). Rather, as electron-phonon coupling is increased, 
the crossover from bound pairs to unbound pairs occu- 
pies a wider range of U values. Comparison with Fig. 
[TJa) shows that total energy at large U typically de- 
creases with increased nearest-neighbor attraction, con- 
sistent with this observation. 

Further evidence for binding of bipolarons can be found 
in the total number of excited phonons shown in Figure 
[2] Here we see that at large Hubbard U the number of 
phonons in the system does not change with respect to 
U, although it is non-zero, because even unbound po- 
larons have phonons associated with them. This can be 
seen in panel (a) the Holstein A increases from 

A = (where there are no phonons) to higher phonon 
coupling where there are a significant number of residual 
phonons at high U. The curve showing the total number 
of phonons also levels out at large U when there is signif- 
icant inter-site coupling and bipolarons large but bound, 
which occurs because there is vanishing on-site compo- 
nent of the wave function in the presence of sufficient U 
and therefore the system is unchanged as U is varied. 

On lowering U in Fig. [5Ja), we see that the num- 
ber of phonons associated with the bipolaron rapidly in- 
creases as the on-site bipolaron forms. This is due to the 
rapid crossover from unbound or inter-site pairs to on- 
site pairs. As U decreases further, the number of phonons 
tends again to a set value dictated by the electron-phonon 
coupling constant and phonon frequency, which also does 
not depend on U. This occurs at sufficient negative U 
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FIG. 2. Number of phonons associated with a singlet 
Hubbard-Holstein bipolaron (a), bipolarons with nearest- 
neighbor interaction strength of 7 = 0.25 (panel b) and 
7 = 0.5 (c). Again, the range of crossover between on-site 
bipolaron and unbound/inter-site bipolarons increases with 
7- 



where the bipolaron is forced into an on-site configu- 
ration. Plots for 7 = 0.25 and 0.5 are shown in Figs, 
[ljb) and[2jc) for the extended Holstein case. Again, an 
increase in nearest neighbor attractive potential visibly 



smooths out transition from bound pairs at low U and the 
unbound or off-site pairs seen at large U. The crossover 
from bound to unbound is stretched over a larger range 
of U consistent with the similar observation in the to- 
tal energy. With increased nearest-neighbor attraction 
the number of associated phonons decreases less dramat- 
ically as on-site pairs form, presumably because an in- 
creased number of phonons are associated with the inter- 
site pairs. 

It is of particular interest to examine the change in 
singlet bipolaron inverse mass as U and A are varied, as 
this can be related to the BEC transition temperature. 
Figure [3] shows that the mass is near constant at large 
U. Bound on-site pairs have a high mass (low inverse 
mass). The inverse mass rapidly decreases at the point 
of binding. We see that the transition from bipolaron to 
polaron starts at lower U, with increased nearest neigh- 
bor interaction [Figs[3]jb) and|3jc)]. The mass decreases 
more slowly with increased Hubbard U for a higher in- 
teraction in accordance with the slow change in associ- 
ated phonons. With high coupling constant A > 1.1 and 
7 = 0.5 we see that the inverse mass has a maximum 
before decreasing and then leveling off, showing that the 
effective mass has a minimum value at intermediate U. 
This phenomenon is not clearly visible in the total en- 
ergy. The reduction in mass is a version of the superlight 
small bipolaron behavio j 4 * 7 * 8 * 31 ! and is achieved when on- 
site and inter-site interactions become comparable in size 
so that bipolarons can hop by contracting and expanding 
through degenerate on-site to inter-site pairs without en- 
ergy penalty. This behavior is significant because small 
mobile bipolarons could form a Bose-Einstein condensate 
with significant transition temperature. We will discuss 
this possibility in the next section. 

Wavefunctions of individual pairs may not overlap if 
bipolarons are to be well defined, so the bipolaron size 
limits the maximum density of particles in a bipolaronic 
material. Figure [4] plots the inverse average singlet bipo- 
laron size against the Hubbard U. Fig. |4][a) shows the 
inverse bipolaron size for the local Holstein interaction. 
As expected for large U, bipolaron pairs unbind and the 
bipolaron size becomes infinite. With increasing coupling 
constant A the average bipolaron size becomes smaller 
(inverse size plotted) as the attractive phonon mediated 
interactions overcome the repulsive Hubbard U. 

In Figs. |4|b) and|4^c), where inter-site interaction is 
turned on, qualitatively different behavior of the bipo- 
laron size can be seen. For weak electron-phonon cou- 
pling A < 1 we see that inverse bipolaron size tends to 
zero at high Hubbard U. The value of U required for 
unbinding increases with 7. With large inter-site inter- 
action of 7 = 0.5, we see that bipolaron unbinding does 
not occur at high U for coupling constant A > 1.1 (within 
the range investigated), instead tending towards a bipo- 
laron size on the order of a lattice constant. At A = 0.8 
and 7 = 0.5, an inflection in the bipolaron size shows the 
transition of on-site bound pairs through off-site pairing 
before the bipolaron completely unbinds on increasing U. 
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FIG. 3. Inverse mass of the singlet bipolaron with nearest- 
neighbor interaction for interaction strength of 7 = (a), 
7 = 0.25 (b) and 7 = 0.5 (c). For large coupling constant A > 
1.1 and 7 = 0.5, we see that the inverse mass has a maximum 
before decreasing and then leveling off. The reduction in mass 
is a version of the superlight small bipolaron behavior, and 
is achieved when on-site and inter-site interactions become 
comparable. Light mobile bipolarons may form a BEC with 
significant transition temperature. 



This is the precursor of the superlight bipolaron behavior 
found at larger A 



FIG. 4. Singlet bipolaron inverse size with nearest-neighbor 
interaction strength of 7 = (a), 7 = 0.25 (b) and 7 = 0.5 
(c). At A = 0.8 and 7 = 0.5, an inflection in the bipolaron size 
shows the transition of on-site bound pairs through inter-site 
pairing before the bipolaron completely unbinds on increasing 
U. This is the precursor of the superlight bipolaron behavior 
found at larger A. 



We conclude this section by examining the infinite 
Hubbard U case. Figure [5] displays inverse bipolaron 
size against increasing electron-phonon coupling constant 
with nearest neighbor interaction strength 7 = 0, 0.25 
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FIG. 5. Singlet bipolaron inverse size with increasing elec- 
tron phonon coupling for infinite Hubbard U, showing the 
transition from unbound to bound states. Nearest neighbor 
strength 7 = 0, 0.25 and 0.5. Inter-site bound pairs are shown 
seen to exist at large A when nearest neighbor interactions are 
switched on. 



and 0.5, showing the crossover from unbound states to 
bound inter-site pairs. It is seen that for 7 = 0.25 and 
7 = 0.5 the system is never fully unbound above A ~ 1.5 
and A ~ 1 respectively. Both have a sharp decrease of 
bipolaron size that tends to a value equal to one lattice 
spacing. Inter-site bound pairs are therefore shown to ex- 
ist at large U and A when nearest neighbor interactions 
are switched on. 



IV. BOSE— EINSTEIN CONDENSATION 

Evidence from the quantum Monte Carlo simulations 
presented in the previous section shows that 3D bipo- 
larons can be simultaneously small and light in the region 
of the parameter space where perturbation theory breaks 
down. In this section, we examine if bosonic charge car- 
riers of this type could form a Bose-Einstein condensate 
(BEC) with a significant transition temperature. 

The BEC transition temperature can be calculated us- 
ing the expression, 



3MH 2 /n 6 \ 2 /3 



B^ BEC 



(10) 



where rib is the number of bosons per site, to** is the 
effective boson mass and a is the lattice constant (here 
we take a to be 4.2A, consistent with the bismuthates). 
An upper bound on the number of bosons per lattice 
site can be established by utilizing the bipolaron size R, 
ni, /a 3 ~ V-R' 3 , where R' is the effective radius of the 
bipolaron. From this we get the following relation, 
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FIG. 6. BEC transition temperature with nearest-neighbor 
interaction strength of 7 = (panel a), 7 = 0.25 (panel b) and 
7 = 0.5 (panel c). The divergence in Tbbc for the Holstein 
case will be suppressed because the bipolaron is only weakly 
bound. On the other hand, high transition temperatures of 
20-30K are seen with nearest neighbor interaction strength 
7 = 0.5 and medium electron-phonon coupling constant where 
bipolarons are well bound. 



Figure [6] is a plot of transition temperatures with (a) 
Holstein and inter-site interaction strengths (b) 7 = 0.25 
and (c) 7 = 0.5. We only plot positive Hubbard U in 



this section, since negative U would be unphysical. The 
effective radius Bl is taken to be 5R (a bipolaron separa- 
tion is 5 bipolaron radii in each direction) corresponding 
to a bipolaron wavefunction overlap of less than 1%. As 
bipolaron densities increase so that bipolarons overlap, 
interaction corrections are expected to reduce transition 
temperatures. 

Several sharp peaks in Tbec are seen in Fig. [6^a) 
for Holstein electron-phonon coupling constants of A = 
0.8, 1.1 and 1.4 (for weaker A, the peaks are at unphysical 
negative U values). These peaks reach values tempera- 
tures of ~ 50K. The regions of high Tbec ar e unstable 
to small variation in U and bipolarons are only weakly 
bound, so thermal fluctuations will suppress transition 
temperatures by breaking up the bipolaron. For weak 
electron-phonon coupling, bipolarons only bind at un- 
physical negative U values, and no BEC is formed for pos- 
itive U. It is interesting to note that although bipolarons 
are formed through electron-phonon coupling, there is a 
wide region of the parameter space where increasing the 
Hubbard U raises the transition temperature. 

BEC transition temperatures are shown for nearest 
neighbor interaction 7 = 0.25 in Figure [6]^b) . Similar 
to the Holstein case, Tbec has a peak at low Hubbard 
U for intermediate A. The peak width increases with 
electron-phonon coupling, but the maximum in the tran- 
sition temperature decreases slightly. A tail appears at 
large U for coupling constant A = 1.4, where the bipo- 
laron becomes bound between sites and properties de- 
pend only weakly on U. QMC simulations are essen- 
tial here, since the regions of high transition temperature 
(where both the electron-phonon coupling and Hubbard 
U are intermediate) can not be accessed using perturba- 
tive techniques. 

Finally Fig. [6](c) plots BEC transition temperatures 
for nearest-neighbor interaction strength 7 = 0.5. The 
strong inter-site coupling completely eradicates the sharp 
peaks in transition temperature associated with the Hol- 
stein interaction, replacing them with broad continuous 
curves with increased Tbec- Bipolarons with medium to 
high A form stable nearest-neighbor pairs with low effec- 
tive masses, leading to superconducting states that have 
significant transition temperatures over wide range of U. 
Bipolarons formed from large electron-phonon coupling 
have larger effective masses, leading to significantly lower 
condensation temperatures. The most interesting point 
here is that for medium-sized coupling constants, bipo- 
laron effective masses are still small when bipolarons are 
bound into small inter-site pairs, resulting in high con- 
densation temperatures of 90 — 120K. Note that the use 
of R' leads to an approximation on the possible TbeCi 
lower transition temperatures are estimated if the up- 
per bound on the distance between bipolarons is larger 
(before inter-boson interactions need to be taken into 
account), and interactions between bipolarons typically 
lower transition temperatures. Again, it is interesting to 
note that the BEC transition temperature can increase 
as U increases, even though the mechanism for binding 



pairs is phonon mediated. Since bipolarons are bound 
at very large U and A (as shown in Figure [4| no break- 
down of the BEC is seen for very large repulsive Coulomb 
interactions, and the transition temperature remains sig- 
nificant. This is important because many oxide materials 
with large electron-phonon interactions also have large U. 
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5 10 15 20 



Hubbard U / 1 

FIG. 7. Effective Hopping energies with nearest-neighbor in- 
teraction strength of 7 = 0.5. Medium electron-phonon cou- 
pling A here are shown to be mobile with effective hopping 
energies in the region of the bare hopping energy. 

To probe the sensitivity of bipolarons to impurities, we 
calculate the effective hopping for the bipolarons, 



Figure [7] shows the effective bipolaron hopping energy 
when 7 = 0.5. For electron-phonon coupling constants 
A = 1.1 and A = 1.4, bipolarons have an effective hopping 
t c g ss 0.7i and t c g sa 0.4t respectively. As we have pre- 
viously showrP, local impurities with energy A = — 4t c g 
are needed to pin polarons to impurities. Therefore, for 
any reasonable impurity size, bipolarons are mobile. 

V. SUMMARY AND CONCLUSIONS 

We have investigated the formation of bipolarons and 
their subsequent Bose - Einstein condensation on a three 
dimensional cubic lattice. A quantum Monte Carlo 
code was employed to investigate regimes of interme- 
diate electron-phonon coupling and Coulomb repulsion, 
and was validated using analytic calculations in the large 
phonon frequency limit. Away from the regions where 
perturbation theories are valid, the effective mass and 
bipolaron radius are consistent with light small bipo- 
larons. 

A consequence of the 3D lattice is that binding of bipo- 
larons is difficult for weak A. Small and mobile bipolarons 
form for intermediate A and inter-site coupling 7 when 
the energies of on-site and inter-site pairs become simi- 
lar. By analyzing the exact numerical results, we have 
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shown that bipolaron condensation temperatures (lead- 
ing to superconductivity) could be up to 90 — 120 K 
for realistic bipolaron densities. Another consequence 
of 3D is that it is more difficult to bind bipolarons to 
impurities. Impurity energies of around — At e s are re- 
quired to localize particles in 3D. Therefore, the light 
bipolaron states are stable against attractive impurity 
levels with energies of up to ~ t as the effective hopping 
has a similar magnitude to the bare electron hopping 
energy. We conclude that stable bipolaron superconduc- 
tors that are insensitive to changes in Hubbard U could 
form in three-dimensional oxides with inter-site electron- 
phonon interactions of intermediate magnitude (that is, 



electron-phonon interactions with a moderate momen- 
tum dependence). Moreover, with sufficiently large inter- 
site electron-phonon coupling, superconductivity could 
be stable at very large values of U, demonstrating that 
Coulomb repulsion is no barrier to bipolaronic supercon- 
ductivity in 3D. 
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Appendix A: High phonon frequency and the UV 
model 



The bipolaron properties can be analytically approx- 
imated in the high phonon frequency (anti-adiabatic) 
limit (Hu) 3> W) by using the result of the Lang-Firsov 
transformation, since if the phonon frequency is very 
large there are no real phonons. Up to a linear shift 
in energy the resulting U — V model is shown as: 



H — t nn ,cl ia .c n ' a + U c i^ c "T c r4 c r4 

nn' a n 

+ ^ ^ V'nn'^noC-ncr^nia'C-n'o' (Al) 



nn' aa' 



The primed sum over V' in the final part of Equation 
(Al) ignores the self- interaction term. The interaction 



terms in this Hamiltonian for on-site interaction and 
nearest neighbor interactions are U' = U — 2WX and 
V', = 2W\ X o(n,n')/ X o(0,0) respectively. 
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Taking the two particle Schrodinger equation, 
[E-e(fe 1 )-e(fe 2 )]x(fei,fe 2 ) 

= u'Y,x(q,k 1 + k 2 -q) 

-V'Y, e~ lkl L E ^1 + fc 2 - g)e iq i (A2) 



e(fe) = -i'^e~ ifeI = -2t'(cosfe :c + cosfe y +cosfc z ](A3) 



where I = {(±1, 0, 0), (0, ±1, 0), (0, 0, ±1)} assuming that 
the lattice constant a = 1 (this will be assumed through- 
out). To simplify the problem we introduce a set of mo- 



mentum dependent values, A(K), 

A (0fitQ) (K) =^x(g,fci + fc2-<?) (A4) 
i 

A l {K) = Y J X{'iM + ki-ci)e i ' 1 - 1 (A5) 



and then substitute them into the Schrodinger equation 
( A2 ) . Rearranging the resulting equation, we obtain an 



expression for $ in terms of K and q, where K — k\ +k2 



x(ki,k 2 ) = 



E - e(fei) - e(k 2 



(A6) 



*<..«-.>- ^V£-r""' '"' 

Exp anding the expression for xili K — q) from Equation 
( A7 ) into a matrix format leads to the following relation: 





L-x 


L m 


L-y 


Ly 


L-z 


Lz 


L x 


L{) + yr 


L2x 


Lx — y 


L x +y 


Lx — z 


L x + z 


L—x 


L-2x 


Lq + 


L — x — y 


L — x+y 


L — x — z 


L — x + z 


Ly 


Ly — x 


Ly+x 


Lo + pr 


L 2 y 


Ly — z 


Ly + z 


L-y 


L — y — x 


L—y-\- X 


L-2y 


Lo + pr 


L-y — z 


L-y + z 


L z 


L z —x 


L z +x 


L Z — y 


Lz+y 


Lo + pr 


L 2 z 


\ L-, 


L — z — X 


L — z + X 


L — Z — y 


L-z+y 


L-2z 


Lq + yr 



( £/'A (0 ,o,o) \ 
-V'A X 

-V'A-x 

-V'Ay 
-V'A-y 

-V'Az 
\ -V'A-z ) 



= 



(A8) 



L p = L qp (E,K)=J2 



E - e{q) - e(K - g) 



(A9) 



At the r point, K = (0,0,0), the additional symme- 
try simplifies analysis of Eq. (A8), since Eq. (A9) 
can be reduced to the following components, L± x = 
L± y = L± z = Lx, L± 2x = L±2y = L± 2z = L 2 , and 
L±x±y = L± y ± z = L± z ± x = L3. 

To diagonalize the problem, a new basis of states with 
s,p and d wave symmetry is introduced; 



An 



Mo. 0,0) 



v/6 



(A* + A-;, +Ay + A-y + A Z + A_ 



1 
1 
1 

: 71 



(A, -A. 

(Ay - A- 

(A z - A_ 



(A10) 



id 2 



A dl 
1 

12 



(A;, + A-z - Ay - A-y 



(A x + A-x + A y + A-y - 2A Z 



2A- 



where the s and d states are symmetric and p states an- 
tisymmetric on inversion through the origin. 

Applying the new basis of A's to the matrix equation 
leads to a simple block diagonal matrix consisting of two 



s states (An, and A s ), three p and two d states. 



/ L - jjt 


6L1 














"\ 


Li 
















































Ap 























Ap 























A d 





V 




















/ U'A \ 
' -V'A S 
-V'A Pl 
-V'A P2 (411) 
-V'A P3 

V -V'Ad, / 



where, 



L s = 

x d 



E 



<3 

--L 

--L 



E - 2e(g) 

2 cos q x (cos q x + cos q y + cos q z ) 
E - 2e(q) 

1 

" V' 
1 

" V' 



(A12) 



L 2 
L 2 



2i a 



The p and d states are diagonalized in the new basis and 
their energies can be calculated directly from, 



X p 

x„ 



(A13) 



The ground state singlet states can be computed from 
solution of the 2x2 matrix in the top left hand corner 
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ofEq. (lATTl) 



Ll Ls + y7 



U'A 
-V'A» 



0, (A14) 



which can be used to calculate energies of the s state by 
taking the determinant, 



Lc 



1 

if' 



1 

r 7 



QL\ = 0. (A15) 



Rearranging Eq. (A15) 



Lo = 



V'(WLl + L s ) + 1 
U'(V'L S + 1) : 



(A16) 



a binary search can be used to determine values of E for 
various values of U' and V . 



of inter-site pairs in the lattice. At large enough V there 
are no unbound states. 

Even if there is no inter-site repulsion, V'/t' = 0, 
strong negative Hubbard values U'/t' = —7.915 are re- 
quired to bind the bipolaron in co ntra st to ID and 2D 
lattices where you need U'/t' = fJ^H This is due to 
the additional degrees of freedom in the cubic structure, 
where electrons are not confined in any direction. For a 
nearest neighbor attraction of V'/t' = 4 Figure [8] shows 
that binding occurs around U' — (actual crossing value 
V'/t' = 3.875), whereas on the square lattice the pres- 
ence of V'/t' leads to off-site pairing at strong positive 
U' . The lack of confinement in 3D mean that bipolaron 
on-site pairing is not guaranteed even with high coupling 
constants. Applying a small nearest-neighbor potential 
in 2D has a much bigger effect on the binding than in 3D 
due to the confinement. 
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FIG. 8. Ground state bipolaron singlet energy computed for 
a UV model on a cubic lattice. A key result here is that it 
takes a large negative Hubbard U to bind on-site bipolarons 
in contrast to ID and 2D systems. At low V' , the bipolaron 
binds for finite U, but for large V' , the unbinding occurs at 
very large (infinite) U since the inter-site bipolaron is stable. 



The solutions found by this method are shown in 
Fig. [H] Typically, there is a smooth transition from 
bound states (diagonal lines) to unbound states (hori- 
zontal lines). For V = the graph depicts a diagonal 
line with constant gradient representing a bound on-site 
pair that at about U'/t' ~ —8 levels off to a horizontal 
line representing unbound polarons in the lattice. With 
increasing nearest-neighbor potential V' the transition 
from bound to unbound states spans a larger range of 
Hubbard U' (curved line), and is related to the presence 



FIG. 9. Exact bipolaron energy computed for a U — V model 
on a cubic lattice. Here we examine the energy at infinite U' 
to observe if there is a critical V' that ensures binding in the 
3D lattice for s, p and d states. 

To understand the qualitative difference in bipolaron 
behavior as V is changed, we evaluated the energy of 
the pair at infinite U . In Fig. [9] we plot the total energy 
(from Equation (A16l) as a function of nearest neighbor 



potential V, showing the binding of the s state as V 
is increased. We also show both p and d states for com- 
pleteness. For the s state, the binding crossover begins at 
a nearest neighbor interaction strength of V'/t' = 5.875. 
After this point the energy curves sharply to an approx- 
imate E oc V' for high V' . The energies of the p and 
d states do not level off at E = —12t', since they are 
excited states. 



